0) Setup

#libzzs---
library(ggplot2)
library(dplyr)
library(ggpubr)
library(car)
library(emmeans)
library(report)
library(DHARMa)
rm(list=ls())
set.seed(666) #Fix kernell

1) Import data & affect factors

setwd("/Users/martinmartin/Downloads/DATA/Github/Field")
sAll<-read.table("Field_G_A_P_learning.csv", sep="",header=T)
sAll$cat<-as.factor(sAll$cat)
sAll$ID<-as.factor(sAll$ID)
sAll$pck<-as.factor(sAll$pck)
sAll$cond<-as.factor(sAll$cond)
sAll$expe<-as.factor(sAll$expe)
sA1<-sAll %>% filter(cond=="GC"|cond=="PC") %>% mutate(cond="CC")
sA2<-sAll %>% filter(cond!="GC"&cond!="PC")
sAll<-bind_rows(sA1,sA2)
sAll$cond<-as.factor(sAll$cond)
levels(sAll$cond)
## [1] "A-" "A+" "AA" "CC" "G-" "G+" "P-" "P+"
c1<-"A 200\u00b5g/L"
c2<-"A 500\u00b5g/L"
c3<-"Control N°1"
c4<-"Control N°2"
c5<-"G 100\u00b5g/L"
c6<-"G 200\u00b5g/L"
c7<-"P 1mg/L"
c8<-"P 10mg/L"
levels(sAll$cond) <- c(c1,c2,c3,c4,c5,c6,c7,c8)

yo<-sAll %>% count(cond,ID)

sAll$cond<-factor(sAll$cond,levels=c("Control N°1",
                                 "A 200\u00b5g/L",
                                 "A 500\u00b5g/L",
                                 "Control N°2",
                                 "G 100\u00b5g/L",
                                 "G 200\u00b5g/L",
                                 "P 1mg/L","P 10mg/L"))
levels(sAll$cond)
## [1] "Control N°1" "A 200µg/L"   "A 500µg/L"   "Control N°2" "G 100µg/L"  
## [6] "G 200µg/L"   "P 1mg/L"     "P 10mg/L"

2) Classification

srawAe<-sAll %>% 
  group_by(cond,expe,ID,cat) %>% 
  mutate(numbering = row_number()) %>% ungroup()

#Define threshold for being "too much low" to be able to dive
l1<- max(srawAe$pos_y)-max(srawAe$pos_y)/8

#Define poscondon to differentiate individuals that did not respond versus respond
lowl1 = -10
upl1 = +10

scaseAe<-srawAe %>% group_by(cond,ID,cat) %>% 
  mutate(threshold = case_when(numbering==min(numbering) & pos_y>l1 ~ 1,TRUE~0)) %>% ungroup()

#Count thresholds and sum into "rep" variable
sfilterAe<-scaseAe %>% 
  group_by(cond,expe,ID,cat) %>%
  summarise(rep=sum(threshold,na.rm=T),
            dY=sum(-dY,na.rm = T),
            pos_y=mean(pos_y,na.rm=T)) %>%  
  filter(rep==0) %>% ungroup()

#Split positions according to limit
sposAe<-sfilterAe %>% 
  mutate(slope=case_when(dY < lowl1 ~ "down",
                         dY>upl1 ~ "up",
                         dY >= lowl1 & dY <= upl1 ~ "flat"))
s5<-sposAe %>% filter(slope=="up") %>% mutate(rep=1)
s5b<-sposAe %>% filter(slope=="flat") %>% mutate(rep=0)
stotAe<-bind_rows(s5,s5b)
rm(s5,s5b,scaseAe,sposAe,srawAe)

srAe<-stotAe %>% 
  group_by(cond,cat) %>% 
  summarise(rep=100*mean(rep,na.rm=TRUE))

stvAe<-sAll %>% 
  group_by(cond,ID,cat) %>% 
  summarise(dY=sum(-dY,na.rm = T),
            pos_y=mean(pos_y,na.rm=T))

complet<-stvAe %>% group_by(cat) %>% count(cat)
filtered<-stotAe %>% group_by(cat) %>% count(cat)
filteredcond<-stotAe %>% group_by(cond,cat) %>% count(cond,cat)

rm(stvAe,complet,filtered)

3) Count number of individuals –> Supp Table T2

jr<-sAll %>% count(cond,ID)
jcat<-sAll %>% count(cond,ID,cat)
yoae<-jcat %>% group_by(cond,ID) %>% count(ID)
yoae1<-yoae %>% filter(cond==c1)
yoae2<-yoae %>% filter(cond==c2)
yoae3<-yoae %>% filter(cond==c3)
yoae4<-yoae %>% filter(cond==c4)
yoae5<-yoae %>% filter(cond==c5)
yoae6<-yoae %>% filter(cond==c6)
yoae7<-yoae %>% filter(cond==c7)
yoae8<-yoae %>% filter(cond==c8)
aebis<-sfilterAe %>% count(cond,ID)
aebis1<-sfilterAe %>% count(cond,ID,cat)
yoaebis<-aebis1 %>% group_by(cond,ID) %>% count(ID)

aeter<-stotAe %>% count(cond,ID)
aeter1<-stotAe %>% count(cond,ID,cat)
yoaeter<-aeter1 %>% group_by(cond,ID) %>% count(ID)

aeae<-left_join(yoae,yoaebis, by=c("cond"="cond","ID"="ID"))
aeae2<-left_join(yoaebis,yoaeter, by=c("cond"="cond","ID"="ID"))

print(c1)
## [1] "A 200µg/L"
print(length(yoae1$ID))
## [1] 30
aec1<-aeae %>% filter(cond==c1)
print("trial per ID")
## [1] "trial per ID"
sum(aec1$n.x)
## [1] 300
sum(aec1$n.y)
## [1] 234
aec1<-aeae2 %>% filter(cond==c1)
sum(aec1$n.y)
## [1] 234
print(c2)
## [1] "A 500µg/L"
print(length(yoae2$ID))
## [1] 29
aec2<-aeae %>% filter(cond==c2)
print("trial per ID")
## [1] "trial per ID"
sum(aec2$n.x)
## [1] 290
sum(aec2$n.y)
## [1] 225
aec2<-aeae2 %>% filter(cond==c2)
sum(aec2$n.y)
## [1] 224
print(c3)
## [1] "Control N°1"
print(length(yoae3$ID))
## [1] 34
aec3<-aeae %>% filter(cond==c3)
print("trial per ID")
## [1] "trial per ID"
sum(aec3$n.x)
## [1] 340
sum(aec3$n.y,na.rm=T)
## [1] 258
aec3<-aeae2 %>% filter(cond==c3)
sum(aec3$n.y)
## [1] 257
print(c4)
## [1] "Control N°2"
print(length(yoae4$ID))
## [1] 44
aec4<-aeae %>% filter(cond==c4)
print("trial per ID")
## [1] "trial per ID"
sum(aec4$n.x)
## [1] 440
sum(aec4$n.y,na.rm=T)
## [1] 367
aec4<-aeae2 %>% filter(cond==c4)
sum(aec4$n.y)
## [1] 362
print(c5)
## [1] "G 100µg/L"
print(length(yoae5$ID))
## [1] 30
aec5<-aeae %>% filter(cond==c5)
print("trial per ID")
## [1] "trial per ID"
sum(aec5$n.x)
## [1] 300
sum(aec5$n.y,na.rm=T)
## [1] 255
aec5<-aeae2 %>% filter(cond==c5)
sum(aec5$n.y)
## [1] 254
print(c6)
## [1] "G 200µg/L"
print(length(yoae6$ID))
## [1] 30
aec6<-aeae %>% filter(cond==c6)
print("trial per ID")
## [1] "trial per ID"
sum(aec6$n.x)
## [1] 300
sum(aec6$n.y,na.rm=T)
## [1] 216
aec6<-aeae2 %>% filter(cond==c6)
sum(aec6$n.y)
## [1] 214
print(c7)
## [1] "P 1mg/L"
print(length(yoae7$ID))
## [1] 31
aec7<-aeae %>% filter(cond==c7)
print("trial per ID")
## [1] "trial per ID"
sum(aec7$n.x)
## [1] 310
sum(aec7$n.y,na.rm=T)
## [1] 241
aec7<-aeae2 %>% filter(cond==c7)
sum(aec7$n.y)
## [1] 240
print(c8)
## [1] "P 10mg/L"
print(length(yoae8$ID))
## [1] 30
aec8<-aeae %>% filter(cond==c8)
print("trial per ID")
## [1] "trial per ID"
sum(aec8$n.x)
## [1] 300
sum(aec8$n.y,na.rm=T)
## [1] 241
aec8<-aeae2 %>% filter(cond==c8)
sum(aec8$n.y)
## [1] 240

4) Evaluate trial bias –> Supp Table T3

print(c1)
## [1] "A 200µg/L"
R1p<-jcat %>% filter(cond==c1) %>% count(cat)
R1r<-stotAe %>% filter(cond==c1) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
## 
##  Chi-squared test for given probabilities
## 
## data:  final1$nz
## X-squared = 13.697, df = 9, p-value = 0.1335
print(c2)
## [1] "A 500µg/L"
R1p<-jcat %>% filter(cond==c2) %>% count(cat)
R1r<-stotAe %>% filter(cond==c2) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
## 
##  Chi-squared test for given probabilities
## 
## data:  final1$nz
## X-squared = 11.879, df = 9, p-value = 0.2202
print(c3)
## [1] "Control N°1"
R1p<-jcat %>% filter(cond==c3) %>% count(cat)
R1r<-stotAe %>% filter(cond==c3) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
## 
##  Chi-squared test for given probabilities
## 
## data:  final1$nz
## X-squared = 5.3133, df = 9, p-value = 0.8062
print(c4)
## [1] "Control N°2"
R1p<-jcat %>% filter(cond==c4) %>% count(cat)
R1r<-stotAe %>% filter(cond==c4) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
## 
##  Chi-squared test for given probabilities
## 
## data:  final1$nz
## X-squared = 11.231, df = 9, p-value = 0.2602
print(c5)
## [1] "G 100µg/L"
R1p<-jcat %>% filter(cond==c5) %>% count(cat)
R1r<-stotAe %>% filter(cond==c5) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
## 
##  Chi-squared test for given probabilities
## 
## data:  final1$nz
## X-squared = 13.565, df = 9, p-value = 0.1387
print(c6)
## [1] "G 200µg/L"
R1p<-jcat %>% filter(cond==c6) %>% count(cat)
R1r<-stotAe %>% filter(cond==c6) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
## 
##  Chi-squared test for given probabilities
## 
## data:  final1$nz
## X-squared = 7.0233, df = 9, p-value = 0.6347
print(c7)
## [1] "P 1mg/L"
R1p<-jcat %>% filter(cond==c7) %>% count(cat)
R1r<-stotAe %>% filter(cond==c7) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
## 
##  Chi-squared test for given probabilities
## 
## data:  final1$nz
## X-squared = 6.2857, df = 9, p-value = 0.711
print(c8)
## [1] "P 10mg/L"
R1p<-jcat %>% filter(cond==c8) %>% count(cat)
R1r<-stotAe %>% filter(cond==c8) %>% count(cat)
Ra1<-left_join(R1p,R1r,by=(c("cat"="cat")))
final1<-Ra1 %>% mutate(nz=n.x-n.y)
chisq.test(final1$nz)
## 
##  Chi-squared test for given probabilities
## 
## data:  final1$nz
## X-squared = 8, df = 9, p-value = 0.5341

5) Plot all trials –> Figure 2) A), B), C)

sf3 <- stotAe %>%
  group_by(cond,cat) %>% 
  summarise(dY=mean(dY,na.rm=T),
            rep=mean(rep,na.rm=T)) %>% ungroup()
sf3$cond<-as.factor(sf3$cond)


st1<-sf3 %>% filter(cond=="Control N°1")
st2<-sf3 %>% filter(cond=="A 200\u00b5g/L")
st3<-sf3 %>% filter(cond=="A 500\u00b5g/L")
st4<-sf3 %>% filter(cond=="Control N°2")
st5<-sf3 %>% filter(cond=="P 1mg/L")
st6<-sf3 %>% filter(cond=="P 10mg/L")
st7<-sf3 %>% filter(cond=="G 100\u00b5g/L")
st8<-sf3 %>% filter(cond=="G 200\u00b5g/L")


print("ATRAZINE")
## [1] "ATRAZINE"
ggplot(NULL) + 
  geom_smooth(data=st1, aes(y=dY, x=cat, group=1, color="Control N°1"), method="gam",
              formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
  geom_smooth(data=st2, aes(y=dY, x=cat, group=1, color="A 200\u00b5g/L"), method="gam",
              formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
  geom_smooth(data=st3, aes(y=dY, x=cat, group=1, color="A 500\u00b5g/L"), method="gam",
              formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
  geom_point(data=st1, aes(y=dY, x=cat, color="Control N°1"), size=3, shape=15) +
  geom_point(data=st2, aes(y=dY, x=cat, color="A 200\u00b5g/L"), size=3, shape=15) +
  geom_point(data=st3, aes(y=dY, x=cat, color="A 500\u00b5g/L"), size=3, shape=15) +
  theme_classic() +
  labs(y="Vertical Distance (mm)", x="Trial", color="Condition") +
  theme(plot.title = element_text(hjust = 0.5, size=22, face="bold"),
        axis.text=element_text(size=20),
        axis.title=element_text(size=20)) +
  scale_x_discrete(name="Trial", 
                   limits=c("1","2","3","4","5","6","7","8","9","10","Test")) +
  scale_color_manual(values = c("Control N°1"="#54C6CC", "A 200\u00b5g/L"="#FFD966", 
                                "A 500\u00b5g/L"="#F0B33E"))

print("PARACETAMOL")
## [1] "PARACETAMOL"
ggplot(NULL) + 
  geom_smooth(data=st4, aes(y=dY, x=cat, group=1, color="Control N°2"), method="gam",
              formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
  geom_smooth(data=st5, aes(y=dY, x=cat, group=1, color="P 1mg/L"), method="gam",
              formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
  geom_smooth(data=st6, aes(y=dY, x=cat, group=1, color="P 10mg/L"), method="gam",
              formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
  geom_point(data=st4, aes(y=dY, x=cat, color="Control N°2"), size=3, shape=15) +
  geom_point(data=st5, aes(y=dY, x=cat, color="P 1mg/L"), size=3, shape=15) +
  geom_point(data=st6, aes(y=dY, x=cat, color="P 10mg/L"), size=3, shape=15) +
  theme_classic() +
  labs(y="Vertical Distance (mm)", x="Trial", color="Condition") +
  theme(plot.title = element_text(hjust = 0.5, size=22, face="bold"),
        axis.text=element_text(size=20),
        axis.title=element_text(size=20)) +
  scale_x_discrete(name="Trial", 
                   limits=c("1","2","3","4","5","6","7","8","9","10","Test")) +
  scale_color_manual(values = c("Control N°2"="#54C6CC", "P 1mg/L"="#00A1EA", 
                                "P 10mg/L"="#0070C0"))

print("GLYPHOSATE")
## [1] "GLYPHOSATE"
ggplot(NULL) + 
  geom_smooth(data=st4, aes(y=dY, x=cat, group=1, color="Control N°2"), method="gam",
              formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
  geom_smooth(data=st7, aes(y=dY, x=cat, group=1, color="G 100\u00b5g/L"), method="gam",
              formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
  geom_smooth(data=st8, aes(y=dY, x=cat, group=1, color="G 200\u00b5g/L"), method="gam",
              formula = y ~ s(x, bs="ps", fx=FALSE, k=-1)) +
  geom_point(data=st4, aes(y=dY, x=cat, color="Control N°2"), size=3, shape=15) +
  geom_point(data=st7, aes(y=dY, x=cat, color="G 100\u00b5g/L"), size=3, shape=15) +
  geom_point(data=st8, aes(y=dY, x=cat, color="G 200\u00b5g/L"), size=3, shape=15) +
  theme_classic() +
  labs(y="Vertical Distance (mm)", x="Trial", color="Condition") +
  theme(plot.title = element_text(hjust = 0.5, size=22, face="bold"),
        axis.text=element_text(size=20),
        axis.title=element_text(size=20)) +
  scale_x_discrete(name="Trial", 
                   limits=c("1","2","3","4","5","6","7","8","9","10","Test")) +
  scale_color_manual(values = c("Control N°2"="#54C6CC", "G 100\u00b5g/L"="#EB8176", 
                                "G 200\u00b5g/L"="#941100"))

6) Plot trials 1 and Test –> Figure 2 D)

sx<-stotAe %>% 
  filter(cat==1|cat==10) %>% 
  group_by(cond,ID,cat) %>%
  summarise(dY=mean(dY,na.rm = T),
            rep=mean(rep,na.rm=T)) %>% ungroup()
sx %>% group_by(cond,cat) %>% count()
## # A tibble: 16 × 3
## # Groups:   cond, cat [16]
##    cond        cat       n
##    <fct>       <fct> <int>
##  1 Control N°1 1        29
##  2 Control N°1 10       26
##  3 A 200µg/L   1        29
##  4 A 200µg/L   10       26
##  5 A 500µg/L   1        23
##  6 A 500µg/L   10       20
##  7 Control N°2 1        40
##  8 Control N°2 10       31
##  9 G 100µg/L   1        28
## 10 G 100µg/L   10       21
## 11 G 200µg/L   1        24
## 12 G 200µg/L   10       22
## 13 P 1mg/L     1        27
## 14 P 1mg/L     10       25
## 15 P 10mg/L    1        27
## 16 P 10mg/L    10       21
dae<-glm(dY~cat*cond, data=sx)
dx<-sx %>% count(cond,cat)
yv <- predict(dae,type="link",se.fit=TRUE)
ydAE<-data.frame(yv$fit,yv$se.fit,sx$cat,sx$cond)
yd2<-ydAE %>% rename(cond=sx.cond) %>% rename(cat=sx.cat)
inv<-family(dae)$linkinv

ggplot(sx, aes(x=cat,y=dY))+
  geom_pointrange(data=yd2,aes(x=cat,y=inv(yv.fit),
                 ymin=inv(yv.fit-1.96*yv.se.fit),
                 ymax=inv((yv.fit+1.96*yv.se.fit)),
                 group=cond,color=cond),size=2,linewidth=3)+
  geom_point(aes(color=cond))+
  theme_classic() +
  labs(y="Vertical Distance (mm)",x="Trial")+
  theme(plot.title = element_text(hjust = 0.5))+
  theme(plot.title = element_text(size=22,face = "bold"))+
  theme(axis.text=element_text(size=20,color="black"),
        axis.title=element_text(size=20,color="black"))+
  scale_color_manual(values=c("#54C6CC","#FFD966","#F0B33E","#54C6CC","#EB8176","#941100",
                             "#00A1EA","#0070C0"))+
  theme(legend.title = element_text(size=20,color="black"),
        legend.text = element_text(size=20))+
  guides(colour=guide_legend(title="Species"))+
  theme(legend.position = "none",
        strip.text = element_text(size = 18))+
  stat_compare_means(comparisons = list(c("1","10")),
                     aes(label = ..p.signif..),
                     bracket.size = 1,size=5)+
  facet_wrap(~cond,nrow=1)

7) Stats

SR1<-sx %>% filter(cond==c1)
mR1<-lmerTest::lmer(dY~cat+(1|ID),data=SR1)
simR1 <- simulateResiduals(fittedModel = mR1, plot = T)

Anova(mR1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: dY
##      Chisq Df Pr(>Chisq)   
## cat 6.9706  1   0.008286 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mR1, pairwise ~ cat,adjust="tukey")
## $emmeans
##  cat emmean   SE df lower.CL upper.CL
##  1     29.1 2.23 53     24.6     33.5
##  10    20.5 2.37 53     15.7     25.2
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate   SE   df t.ratio p.value
##  cat1 - cat10     8.57 3.26 27.7   2.631  0.0137
## 
## Degrees-of-freedom method: kenward-roger
SR2<-sx %>% filter(cond==c2)
mR2<-lmerTest::lmer(dY~cat+(1|ID),data=SR2)
simR2 <- simulateResiduals(fittedModel = mR2, plot = T)

Anova(mR2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: dY
##      Chisq Df Pr(>Chisq)  
## cat 4.0355  1    0.04455 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mR2, pairwise ~ cat,adjust="tukey")
## $emmeans
##  cat emmean   SE   df lower.CL upper.CL
##  1     31.5 2.67 40.9     26.1     36.9
##  10    24.0 2.88 41.0     18.2     29.8
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate   SE   df t.ratio p.value
##  cat1 - cat10     7.46 3.78 23.2   1.974  0.0604
## 
## Degrees-of-freedom method: kenward-roger
SR3<-sx %>% filter(cond==c3)
mR3<-lmerTest::lmer(dY~cat+(1|ID),data=SR3)
simR3 <- simulateResiduals(fittedModel = mR3, plot = T)

Anova(mR3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: dY
##      Chisq Df Pr(>Chisq)    
## cat 11.613  1  0.0006549 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mR3, pairwise ~ cat,adjust="tukey")
## $emmeans
##  cat emmean   SE   df lower.CL upper.CL
##  1     29.7 2.03 49.9     25.6     33.8
##  10    21.4 2.14 51.3     17.1     25.7
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate   SE   df t.ratio p.value
##  cat1 - cat10     8.28 2.46 26.7   3.373  0.0023
## 
## Degrees-of-freedom method: kenward-roger
SR4<-sx %>% filter(cond==c4)
mR4<-lmerTest::lmer(dY~cat+(1|ID),data=SR4)
simR4 <- simulateResiduals(fittedModel = mR4, plot = T)

Anova(mR4)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: dY
##      Chisq Df Pr(>Chisq)   
## cat 7.9536  1   0.004799 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mR4, pairwise ~ cat,adjust="tukey")
## $emmeans
##  cat emmean   SE df lower.CL upper.CL
##  1     26.5 1.85 69     22.9     30.2
##  10    18.7 2.12 69     14.4     22.9
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate   SE   df t.ratio p.value
##  cat1 - cat10     7.88 2.81 38.5   2.801  0.0079
## 
## Degrees-of-freedom method: kenward-roger
SR5<-sx %>% filter(cond==c5)
mR5<-lmerTest::lmer(dY~cat+(1|ID),data=SR5)
simR5 <- simulateResiduals(fittedModel = mR5, plot = T)

Anova(mR5)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: dY
##      Chisq Df Pr(>Chisq)
## cat 1.1381  1     0.2861
emmeans(mR5, pairwise ~ cat,adjust="tukey")
## $emmeans
##  cat emmean   SE df lower.CL upper.CL
##  1     22.7 2.54 47     17.6     27.8
##  10    18.5 2.97 47     12.6     24.5
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate   SE   df t.ratio p.value
##  cat1 - cat10     4.14 3.91 25.2   1.059  0.2995
## 
## Degrees-of-freedom method: kenward-roger
SR6<-sx %>% filter(cond==c5)
mR6<-lmerTest::lmer(dY~cat+(1|ID),data=SR6)
simR6 <- simulateResiduals(fittedModel = mR6, plot = T)
Anova(mR6)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: dY
##      Chisq Df Pr(>Chisq)
## cat 1.1381  1     0.2861
emmeans(mR6, pairwise ~ cat,adjust="tukey")
## $emmeans
##  cat emmean   SE df lower.CL upper.CL
##  1     22.7 2.54 47     17.6     27.8
##  10    18.5 2.97 47     12.6     24.5
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate   SE   df t.ratio p.value
##  cat1 - cat10     4.14 3.91 25.2   1.059  0.2995
## 
## Degrees-of-freedom method: kenward-roger
SR7<-sx %>% filter(cond==c7)
mR7<-lmerTest::lmer(dY~cat+(1|ID),data=SR7)
simR7 <- simulateResiduals(fittedModel = mR7, plot = T)

Anova(mR7)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: dY
##      Chisq Df Pr(>Chisq)   
## cat 6.9368  1   0.008444 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(mR7, pairwise ~ cat,adjust="tukey")
## $emmeans
##  cat emmean   SE   df lower.CL upper.CL
##  1     25.0 2.40 45.7     20.1     29.8
##  10    17.6 2.49 46.9     12.6     22.6
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate   SE   df t.ratio p.value
##  cat1 - cat10     7.37 2.81 25.1   2.618  0.0148
## 
## Degrees-of-freedom method: kenward-roger
SR8<-sx %>% filter(cond==c8)
mR8<-lmerTest::lmer(dY~cat+(1|ID),data=SR8)
simR8 <- simulateResiduals(fittedModel = mR8, plot = T)

Anova(mR8)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: dY
##      Chisq Df Pr(>Chisq)
## cat 2.0166  1     0.1556
emmeans(mR8, pairwise ~ cat,adjust="tukey")
## $emmeans
##  cat emmean   SE   df lower.CL upper.CL
##  1     26.4 2.51 37.6     21.3     31.4
##  10    22.6 2.76 42.7     17.1     28.2
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate   SE   df t.ratio p.value
##  cat1 - cat10     3.73 2.65 21.8   1.408  0.1732
## 
## Degrees-of-freedom method: kenward-roger